England

# remotes::install_github("mrc-ide/COVIDCurve", ref = "develop")
library(COVIDCurve)
library(tidyverse)
source("R/extra_plotting_functions.R")
source("R/my_themes.R")

Results

Age Estimates

#...................... 
# read in data
#......................
modout <- readRDS("results/ModFits/GBR_age_rung50_burn10000_smpl10000.RDS")

Log Plot

modout$mcmcout$output %>%
  dplyr::filter(rung == "rung1") %>%
  dplyr::filter(stage == "sampling") %>%
  dplyr::select(c("chain", "iteration", "loglikelihood", "logprior")) %>%
  tidyr::gather(., key = "like", value = "val", 3:4) %>%
  ggplot() +
  geom_line(aes(x = iteration, y = val, color = like), size = 0.1, alpha = 0.8) +
  facet_wrap(chain~like, scales = "free_y") +
  theme_bw()

Quick MCMC Diagnostics

quickout <- quick_mc_diagnostics(modout)
quickout

IFRs & Incidence Curve

agekey <- modout$inputs$IFRmodel$IFRdictkey %>% 
  dplyr::rename(param = Strata)
#...................... 
# get ifrs 
#......................
ifrs <- COVIDCurve::get_cred_intervals(IFRmodel_inf = modout, whichrung = paste0("rung", 1),
                                       what = "IFRparams", by_chain = FALSE)
ifrs <- dplyr::left_join(agekey, ifrs, by = "param") %>% 
  dplyr::mutate(param = factor(param, levels = paste0("ma", 1:5))) %>% 
  dplyr::arrange()

#...................... 
# get crude ifrs 
#......................
dscdat <- readRDS("data/derived/descriptive_results_datamap.RDS") 
GBRcrudedat <- dscdat %>% 
  dplyr::filter(study_id == "GBR2" & breakdown == "ageband")  %>%
  dplyr::select(c("study_id", "plotdat")) %>%
  tidyr::unnest(cols = "plotdat")
ageplotdat <- GBRcrudedat %>%
  dplyr::filter(seromidpt == obsday) %>%
  dplyr::mutate(infxns = popn * seroprev,
                crudeIFR =  cumdeaths/infxns) %>% 
  dplyr::select(c("ageband", "crudeIFR")) 


#...................... 
# get incidence curve
#......................
infxncurve <- COVIDCurve::draw_posterior_infxn_cubic_splines(IFRmodel_inf = modout,
                                                             dwnsmpl = 1e3,
                                                             by_chain = FALSE, 
                                                             by_strata = TRUE)
#...................... 
# make ifr and incidence plots
#......................
plot1 <- ggplot() +
  geom_pointrange(data = ifrs, aes(x = ageband, ymin = LCI, ymax = UCI, y = median), 
                  color = "#969696", size = 1.2) +
  geom_point(data = ageplotdat, aes(x = ageband, y = crudeIFR), color = "#000000", size = 2) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, vjust = 0.90, hjust= 1, face = "bold"),
        legend.position = "right") +
  xlab("") + ylab("Median (95% CIs)") +
  labs(caption = "Grey and black points represent model fits and crude IFR estimates, respectively")
plot2 <- infxncurve$plotObj
# main plot
main_plotObj <- cowplot::plot_grid(plot1, plot2, ncol = 1, nrow = 2)

# quick out
jpgsnapshot("figures/GBR2_ifr_plot_age.jpg", plot = main_plotObj)
plot(main_plotObj)

ifrs %>% 
  dplyr::mutate_if(is.numeric, round, 2) %>% 
  pretty_DT_tab(.)
Overall IFR
# population denom
demog <- modout$inputs$IFRmodel$demog %>% 
  dplyr::rename(param = Strata)
demog <- demog %>% 
  dplyr::left_join(., agekey)

gIFR <- sum(ifrs$median * demog$popN/sum(demog$popN))

tibble::tibble("Overall IFR" = sprintf("%1.2f%%", 100*gIFR)) %>% 
  knitr::kable(.)
Overall IFR
3.27%

Seroprevalence Observed vs. Inferred True Prev

seropnts <- COVIDCurve::draw_posterior_sero_curves(IFRmodel_inf = modout,
                                                   dwnsmpl = 1e3,
                                                   by_chain = F)

serocurvedat <- seropnts %>% 
  dplyr::select(c("sim", "ObsDay", dplyr::starts_with("RG_"), dplyr::starts_with("inf_"))) %>% 
  tidyr::gather(., key = "seroprev_strata_lvl", value = "serocount", 3:ncol(.)) %>% 
  dplyr::mutate(seroprevlvl = ifelse(stringr::str_detect(seroprev_strata_lvl, "RG_"), "RG", "IF"),
                param = stringr::str_extract(seroprev_strata_lvl, "ma[0-9]+")) %>% 
  dplyr::group_by_at(c("sim", "ObsDay", "seroprevlvl", "param")) %>% 
  dplyr::summarise(serocount = sum(serocount, na.rm = T) ) %>% 
  dplyr::left_join(., agekey, by = "param") %>% 
  dplyr::left_join(., demog, by = "ageband") %>% 
  dplyr::mutate(seroprev = serocount/popN) 

# crude data
SeroPrevObs <- GBRcrudedat %>% 
  dplyr::select(c("obsdaymin", "obsdaymax", "ageband", "seroprev")) %>% 
  dplyr::filter(!duplicated(.)) %>% 
  dplyr::left_join(., demog, by = "ageband") %>% 
  dplyr::mutate(obsmidday = (obsdaymin + obsdaymax)/2)
serocurvedat %>% 
  ggplot() +
  geom_line(aes(x = ObsDay, y = seroprev, color = seroprevlvl), alpha = 0.5) +
  geom_rect(data = SeroPrevObs, aes(xmin = obsdaymin, xmax = obsdaymax, ymin = -Inf, ymax = Inf), 
            color = "#d9d9d9", fill = "#d9d9d9", alpha = 0.8) +
  geom_point(data = SeroPrevObs, aes(x = obsmidday, y = seroprev, group = ageband),
             color = "#000000", size = 1.2) +
  facet_wrap(.~ageband) +
  scale_color_manual("Seroprev. \n Adjustment", values = c("#FFD301", "#246BCF"),
                     labels = c("Inferred 'Truth'", "Inferred 'Observed' - \n Rogan-Gladen Corrected")) +
  labs(caption = "Grey box is observed seroprevalence across study period") +
  xyaxis_plot_theme + 
  theme(axis.text.x = element_text(angle = 45, vjust = 0.90, hjust= 1, face = "bold"))

Posterior Predictive Check

#...................... 
# get posterior draws
#......................
postdat <- COVIDCurve::posterior_check_infxns_to_death(IFRmodel_inf = modout,
                                                       dwnsmpl = 1e3,
                                                       by_chain = FALSE)
postdat_long <- postdat %>%
  dplyr::select(c("sim", "time", dplyr::starts_with("deaths"))) %>%
  tidyr::gather(., key = "param", value = "deaths", 3:ncol(.)) %>%
  dplyr::mutate(param = gsub("deaths_", "", param)) %>% 
  dplyr::left_join(., y = agekey, by = "param") 

#......................
# get deaths posterior pred check
#......................
datclean <-  modout$inputs$IFRmodel$data$obs_deaths %>%
  dplyr::filter(Deaths != -1) %>% 
  dplyr::left_join(., modout$inputs$IFRmodel$IFRdictkey, by = "Strata")



#...................... 
# make plot
#......................
ppc_plot_age <- ggplot() +
  geom_line(data = postdat_long, aes(x= time, y = deaths, group = ageband), 
            size = 1.2, color = "#bdbdbd") +
  geom_line(data = datclean, aes(x=ObsDay, y = Deaths, group = ageband), 
            color = "#3182bd") +
  facet_wrap(.~ageband) +
  theme_bw() +
  ggtitle("Posterior Predictive Check", subtitle = "Grey Lines are Draws from Posterior, Blue Line is Observed Data (from ECDC)") +
  theme(plot.title = element_text(hjust = 0.5, vjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5, vjust = 0.5))

jpgsnapshot("figures/GBR2_ppc_plot_age.jpg", plot = ppc_plot_age)
plot(ppc_plot_age)

Regional Estimates

#...................... 
# read in data
#......................
modout <- readRDS("results/ModFits/GBR_rgn_rung50_burn10000_smpl10000.RDS")

Log Plot

modout$mcmcout$output %>%
  dplyr::filter(rung == "rung1") %>%
  dplyr::filter(stage == "sampling") %>%
  dplyr::select(c("chain", "iteration", "loglikelihood", "logprior")) %>%
  tidyr::gather(., key = "like", value = "val", 3:4) %>%
  ggplot() +
  geom_line(aes(x = iteration, y = val, color = like), size = 0.1, alpha = 0.8) +
  facet_wrap(chain~like, scales = "free_y") +
  theme_bw()

Quick MCMC Diagnostics

quickout <- quick_mc_diagnostics(modout)
quickout

IFRs & Incidence Curve

rgnkey <- modout$inputs$IFRmodel$IFRdictkey %>% 
  dplyr::rename(param = Strata)
#...................... 
# get ifrs 
#......................
ifrs <- COVIDCurve::get_cred_intervals(IFRmodel_inf = modout, whichrung = paste0("rung", 1),
                                       what = "IFRparams", by_chain = FALSE)
ifrs <- dplyr::left_join(rgnkey, ifrs, by = "param") %>% 
  dplyr::mutate(param = factor(param, levels = paste0("ma", 1:3))) %>% 
  dplyr::arrange()

#...................... 
# get crude ifrs 
#......................
dscdat <- readRDS("data/derived/descriptive_results_datamap.RDS") 
GBRcrudedat <- dscdat %>% 
  dplyr::filter(study_id == "GBR2" & breakdown == "region")  %>%
  dplyr::select(c("study_id", "plotdat")) %>%
  tidyr::unnest(cols = "plotdat")
rgnplotdat <- GBRcrudedat %>%
  dplyr::filter(seromidpt == obsday) %>%
  dplyr::mutate(infxns = popn * seroprev,
                crudeIFR =  cumdeaths/infxns) %>% 
  dplyr::select(c("region", "crudeIFR")) 


#...................... 
# get incidence curve
#......................
infxncurve <- COVIDCurve::draw_posterior_infxn_cubic_splines(IFRmodel_inf = modout,
                                                             dwnsmpl = 1e3,
                                                             by_chain = FALSE, 
                                                             by_strata = TRUE)
#...................... 
# make ifr and incidence plots
#......................
plot1 <- ggplot() +
  geom_pointrange(data = ifrs, aes(x = region, ymin = LCI, ymax = UCI, y = median), 
                  color = "#969696", size = 1.2) +
  geom_point(data = rgnplotdat, aes(x = region, y = crudeIFR), color = "#000000", size = 2) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, vjust = 0.90, hjust= 1, face = "bold"),
        legend.position = "right") +
  xlab("") + ylab("Median (95% CIs)") +
  labs(caption = "Grey and black points represent model fits and crude IFR estimates, respectively")
plot2 <- infxncurve$plotObj
# main plot
main_plotObj <- cowplot::plot_grid(plot1, plot2, ncol = 1, nrow = 2)

# quick out
jpgsnapshot("figures/GBR2_ifr_plot_rgn.jpg", plot = main_plotObj)
plot(main_plotObj)

ifrs %>% 
  dplyr::mutate_if(is.numeric, round, 2) %>% 
  pretty_DT_tab(.)
Overall IFR
# population denom
demog <- modout$inputs$IFRmodel$demog %>% 
  dplyr::rename(param = Strata)
demog <- demog %>% 
  dplyr::left_join(., rgnkey)

gIFR <- sum(ifrs$median * demog$popN/sum(demog$popN))

tibble::tibble("Overall IFR" = sprintf("%1.2f%%", 100*gIFR)) %>% 
  knitr::kable(.)
Overall IFR
1.28%

Seroprevalence Observed vs. Inferred True Prev

seropnts <- COVIDCurve::draw_posterior_sero_curves(IFRmodel_inf = modout,
                                                   dwnsmpl = 1e3,
                                                   by_chain = F)

serocurvedat <- seropnts %>% 
  dplyr::select(c("sim", "ObsDay", dplyr::starts_with("RG_"), dplyr::starts_with("inf_"))) %>% 
  tidyr::gather(., key = "seroprev_strata_lvl", value = "serocount", 3:ncol(.)) %>% 
  dplyr::mutate(seroprevlvl = ifelse(stringr::str_detect(seroprev_strata_lvl, "RG_"), "RG", "IF"),
                param = stringr::str_extract(seroprev_strata_lvl, "ma[0-9]+")) %>% 
  dplyr::group_by_at(c("sim", "ObsDay", "seroprevlvl", "param")) %>% 
  dplyr::summarise(serocount = sum(serocount, na.rm = T) ) %>% 
  dplyr::left_join(., rgnkey, by = "param") %>% 
  dplyr::left_join(., demog, by = "region") %>% 
  dplyr::mutate(seroprev = serocount/popN) 

# crude data
SeroPrevObs <- GBRcrudedat %>% 
  dplyr::select(c("obsdaymin", "obsdaymax", "region", "seroprev")) %>% 
  dplyr::filter(!duplicated(.)) %>% 
  dplyr::left_join(., demog, by = "region") %>% 
  dplyr::mutate(obsmidday = (obsdaymin + obsdaymax)/2)
serocurvedat %>% 
  ggplot() +
  geom_line(aes(x = ObsDay, y = seroprev, color = seroprevlvl), alpha = 0.5) +
  geom_rect(data = SeroPrevObs, aes(xmin = obsdaymin, xmax = obsdaymax, ymin = -Inf, ymax = Inf), 
            color = "#d9d9d9", fill = "#d9d9d9", alpha = 0.8) +
  geom_point(data = SeroPrevObs, aes(x = obsmidday, y = seroprev, group = region),
             color = "#000000", size = 1.2) +
  facet_wrap(.~region) +
  scale_color_manual("Seroprev. \n Adjustment", values = c("#FFD301", "#246BCF"),
                     labels = c("Inferred 'Truth'", "Inferred 'Observed' - \n Rogan-Gladen Corrected")) +
  labs(caption = "Grey box is observed seroprevalence across study period") +
  xyaxis_plot_theme + 
  theme(axis.text.x = element_text(angle = 45, vjust = 0.90, hjust= 1, face = "bold"))

Posterior Predictive Check

#...................... 
# get posterior draws
#......................
postdat <- COVIDCurve::posterior_check_infxns_to_death(IFRmodel_inf = modout,
                                                       dwnsmpl = 1e3,
                                                       by_chain = FALSE)
postdat_long <- postdat %>%
  dplyr::select(c("sim", "time", dplyr::starts_with("deaths"))) %>%
  tidyr::gather(., key = "param", value = "deaths", 3:ncol(.)) %>%
  dplyr::mutate(param = gsub("deaths_", "", param)) %>% 
  dplyr::left_join(., y = rgnkey, by = "param") 

#......................
# get deaths posterior pred check
#......................
datclean <-  modout$inputs$IFRmodel$data$obs_deaths %>%
  dplyr::filter(Deaths != -1) %>% 
  dplyr::left_join(., modout$inputs$IFRmodel$IFRdictkey, by = "Strata")



#...................... 
# make plot
#......................
ppc_plot_rgn <- ggplot() +
  geom_line(data = postdat_long, aes(x= time, y = deaths, group = region), 
            size = 1.2, color = "#bdbdbd") +
  geom_line(data = datclean, aes(x=ObsDay, y = Deaths, group = region), 
            color = "#3182bd") +
  facet_wrap(.~region) +
  theme_bw() +
  ggtitle("Posterior Predictive Check", subtitle = "Grey Lines are Draws from Posterior, Blue Line is Observed Data (from ECDC)") +
  theme(plot.title = element_text(hjust = 0.5, vjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5, vjust = 0.5))

jpgsnapshot("figures/GBR2_ppc_plot_rgn.jpg", plot = ppc_plot_rgn)
plot(ppc_plot_rgn)